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Abstract 

Improvements beyond the primitive approximation in the path integral Monte Carlo method are 
explored both in a model problem and in real systems. Two different strategies are studied: the 
Richardson extrapolation on top of the path integral Monte Carlo data and the Takahashi-Imada 
action. The Richardson extrapolation, mainly combined with the primitive action, always reduces 
the number-of-beads dependence, helps in determining the approach to the dominant power law 
behavior, and all without additional computational cost. The Takahashi-Imada action has been 
tested in two hard-core interacting quantum liquids at low temperature. The results obtained show 
that the fourth-order behavior near the asymptote is conserved, and that the use of this improved 
action reduces the computing time with respect to the primitive approximation. 

PACS numbers: 31.15.Kb,02.70.Ss 
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I. INTRODUCTION 



Path integral Monte Carlo (PIMC) has become in the last decade a standard tool for 
studying quantum liquids and solids at finite temperature. 1-3 The PIMC method allows for 
the calculation of the quantum-statistical partition function, and from it, thermodynamic 
functions like the internal energy or the specific heat. The basis on which PIMC rests 
is the analytical continuation to imaginary time of Feynman's path-integral formalism of 
Quantum Mechanics. 4-6 The resulting convolutive property for the thermal density matrix 
makes accessible the quantum physics of low temperatures starting from density matrices 
at higher temperatures. At sufficiently high temperature, commutator terms between the 
kinetic and potential operators defining the Hamiltonian of the system can be disregarded, 
and then, the density matrix factorizes as a product of the kinetic and potential parts. This 
is known as the primitive approximation (PA) and the unbiased convergence to a desired 
lower temperature by increasing the number M of convolutive terms is guaranteed by the 
Trotter formula. 7 

The expectation value of an operator through the thermal density matrix is written in 
this language as a multidimensional integral appropriate for a Monte Carlo interpretation. 
The quantum problem can be mapped to a classical one of interacting polymers, each one 
representing an atom. 8,9 In the action, beads of different atoms but with the same index 
interact through the potential, and all the beads of the same atom form a closed chain of 
springs. At the level of the PA and avoiding the specific quantum statistics of the atoms or 
molecules, relevant at temperatures near the absolute zero, the PIMC method looks rather 
simple. However, in the implementation of the algorithm there are technical aspects which 
require of some additional effort. The main point is that, as the number of beads of the 
polymeric chains increases rapidly with decreasing temperature, the sampling of the proba- 
bility distribution function becomes almost impossible by means of individual movements. 
Several mechanisms for performing smart collective movements of the beads have been pro- 
posed and proved to be able of eliminating the slowing down present in the single-bead 
schemes. Worth mentioning are the staging 10-12 and bisection 2 ' 13 which exploit the fact that 
the free action (kinetic part in an interacting system) can be sampled by non-rejection meth- 
ods, and therefore, only the potential part requires a Metropolis step. In a large number of 
PIMC simulations, these collective sampling schemes are enough to reach a reasonable de- 
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gree of efficiency. Nevertheless, calculations of quantum fluids at extremely low temperature 
require chains composed by several thousands of beads, feature which makes the problem 
more involved. 

The number of beads required to achieve the asymptotic regime, at a given tempera- 
ture, depends on the approximation taken for the action. It has been proved that PA is 
accurate to order (/3/M) 2 , with f3 the inverse temperature and M the number of beads. 14 ' 15 
Focusing the analysis on PIMC, and thus not considering the specific developments for the 
Fourier PIMC method, 16 the improvement of the action can be pursued following at least 
two different alternatives. In the first one, the PA is substituted by the pair-product action 
in which the basic step on the chain is constructed as the exact action for two isolated 
atoms. 10 This method is more efficient than the PA, mainly when dealing with hard-sphere- 
like systems, and it is the only method up to date which has proved to be accurate enough 
to reach the superfluid regime of liquid 4 He. 2 The second possibility relies on the inclusion 
in the action of the lowest-order corrections to the primitive approximation for the expo- 
nential of the Hamiltonian e~^ H . To this end, Takahashi and Imada 17 and later on, and 
independently, Li and Broughton 18 found a manageable expression for the trace accurate 
to order (/3/M) 4 . In this approximation, hereafter referred as TIA, the beads of different 
atoms interact through an effective potential which is the sum of the interatomic potential 
and a /9-dependent term containing the double commutator [[V, K],V], with V and K the 
potential and kinetic operators, respectively. Recently, Jang et al. 19 have used in PIMC a 
genuine fourth-order factorization according to the proposal of Suzuki 20 and the posterior 
developments of Chin. 21 ' 22 However, the accuracy achieved in several test problems by Jang 
et al. 19 is comparable to the one reached using the simpler TIA. 

In the present paper we analyze the f3/M dependence of the energy using both PA and 
TIA for a test model and for two strongly interacting fluids at low temperature, i.e., Ne and 
4 He. In the analysis of the PA results, we introduce the Richardson extrapolation. 23 The 
efficiency of this numerical scheme for extrapolation is well known in fields like integration 
or solution of differential equations. We show that its use in PIMC can also be useful in 
two directions. First, its simple use helps in determining when the expected law ((3/M) 2 
is reached; second, the successive extrapolations to M — > oo follow a nearly fourth-order 
dependence and then they always approach faster to the asymptote than the PA. The second 
main objective of the work is the study of the TIA when applied to hard-core interacting 
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systems in the regime of low temperatures since the applications of TIA to these situations 
is rather scarce. 24-27 Previous tests of the TIA to systems not so dense and not so cold have 
shown the expected fourth-order accuracy. 24 In fact, the achievement of a fourth-order law 
for hard-core potentials using TIA has been some times questioned 2 but never thoroughly 
analyzed. Our present results show that in the regime studied the accuracy obtained is the 
expected one and the global efficiency of TIA is larger than the one of the standard PA. 
On the other hand, the use of the Richardson extrapolation on top of the TIA results has 
proved not to be so helpful, at least in the temperatures where the two real systems have 
been studied. 

The rest of the paper is organized as follows: Section II contains the formalism of the PA 
and the TIA, including the expressions for the thermodynamic and centroid-virial estimators 
for the energy. The application of the Richardson extrapolation to the PIMC calculations 
is also discussed. In Sec. Ill, results for the three systems under study are presented with 
special emphasis to the {(3/M) dependence of the energy obtained using the two models for 
the action. Finally, Sec. IV comprises a brief summary and the main conclusions of the 
work. 



II. FORMALISM 

The partition function of a system described by a Hamiltonian H at a temperature 
T = 1/(3 is 

Z = Tr e~ pS . (1) 
We consider a system composed by N particles with H = K + V , the kinetic operator being 

fi 2 N 

i=i 

and the potential one 

N 

V > = E^)' ( 3 ) 

i<j 

where pairwise interactions are assumed. The PA is the lowest order term in powers of 
((3/M) of the exponential of the sum of the two operators K and V, 

M 



Tr e -/3(if+t>) ~ Tr ^ e -(^/M) e -(pv/M) j _ ^ 
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In fact, in the limit where M — > oo the previous expression becomes exact, yielding the 
Trotter formula. With the decomposition (4), and restricting the analysis to distinguishable 
particles, the partition function becomes 

/M 
dH\ . . . dR M JJp PA ( R a , R a +i) , 
(5) 

a=l 

with R = {t*i, . . . , r N } and Rm+i = Ri- Greek and Latin indexes are used throughout the 
paper to denote beads and particles, respectively. The primitive factorization (4), which 
decouples the kinetic and potential operators, leads to the well-known PA, 

/ Mm \3V/2 f N Mm p N 

p PA (R a ,R a+1 ) = (j^aj exp|-J]^(r a)i -r Q+M ) 2 - j^Yl V ^ a ^j " (6) 

Takahashi and Imada 17 proved that it is possible to extend the accuracy of the action up 
to fourth order for the trace by substituting in Eq. (4) the operator V by another one W 
given by 

i<j v 7 i=l 

The second term on the r.h.s of Eq. (7) corresponds to the double commutator [[V, K], V] 
appearing in the development of the exponential e~@ H . The classical-like force is defined 

as 

N 

Therefore, the partition function using the TIA is written like the PA one (5) just substi- 
tuting the interatomic potential V(r a jj) in Eq. (6) by W(r a jj), 



< n n N (Mm\ 3N/2 f *Mm f 



(9) 



\i<j x ' i=l 

First derivatives of the partition function allow for the estimation of the total and partial 
energies of the iV-body system, 5 

E 1 dZ 

N ~ ~Wz~d(3 ( ' 

N N(3Z dm v ; 
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and V/N = E/N — K/N . The potential energy, and any other coordinate operators, can 
also be obtained through the general relation 17 

1 1 dZ(V + AO) 



O(R) 



(12) 

A=0 



PZ(V) dX 

Equation (11) leads to the thermodynamic estimators for the kinetic energy. At the level of 
PA, 



Ki\ 3M " " Mm 



N 28 ^ ^ 2N(3 2 h 2 

^ a=l i=l ' 

In the case of the TIA the effective potential W depends on the mass of the particles and 
the temperature. Therefore, explicit terms depending on the interatomic potential appear 
in the estimation of the kinetic energy. Explicitly, 

Kf 1A _ , 1 h 2 (3 2 _ w 

3 

a=l i=l 



The potential energy in the PA case corresponds to a direct estimation of V(r), 

1 M N 



N NM 

a=l i<j 



whereas the application of Eq. (12) for the TIA produces the same additional force term of 
Eq. (14), multiplied by a factor of two, 



.2 n 2 m N 



VTIA _ VPA , 1 h P Y^^lp |2 / lfi \ 

iV AT + 12 m NM^ ^ l aM ' 1 j 

The unbiased convergence of both the PA and TIA factorizations to the exact energy when 
M increases is granted by the Trotter formula. 7 However, the thermodynamic estimator for 
both approaches (13,14) presents the drawback of a statistical variance which increases with 
the number of beads M. 28,29 To overcome this problem one can derive another estimator 
for the energy by exploiting the invariance of the partition function under a rescaling of the 
particle positions (r* — > Ar*). This leads to the so called virial estimator 28 in which the 
kinetic energy is obtained from derivatives of the interatomic potential V(r). It has been 
proved that, contrary to the thermodynamic estimator, the variance of the virial estimator 
is roughly constant as a function of M. 28,29 We have chosen in the present work the centroid 
version of the virial estimator. In the PA, 

1 1 M N 



N 28 ' 2 NM 

( a=l i=l 
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with r 0ji = ( r^j + • • • + r a> i)/M the centroid position of atom i. By applying the same 
method to the TIA, two new terms appear requiring up to the second derivative of V(r), 

a=l i=l L j^j 

(18) 

where {a, 6} stand for the Cartesian coordinates and an implicit summation over repeated 
indices has been assumed. The tensor T is explicitly written as 



T(a,i,j) b a = 



dr n a ri. „■„■ dr 2 



r <*,ij r a,ij 



+ —2 " 1Z2~^- ( l9 ) 



a>*J ' a,ij a,ij 

The PIMC method is the most efficient theoretical tool to deal with a microscopic de- 
scription of systems at low temperatures where quantum effects are unavoidable. With the 
only external assumption of the interatomic potentials, PIMC provides exact results at a 
given temperature with a self-contained adjustment of the number of terms (beads M) in the 
action to the quanticity of the system: from M — 1, which corresponds to classical Monte 
Carlo (each particle is a point in the configuration space), to values of M (each particle is 
represented by a closed chain of beads) larger enough to fulfill the Trotter formula. 7 It is 
hence essential to remove from the calculation the bias coming from the use of a discrete 
value for M, or at least to reduce it to the level of the typical statistical error. In order 
to help in this analysis we have used the well-known Richardson extrapolation; 23 its use in 
numerical integration, derivation or differential equations permits to achieve high accuracy 
by using low-order formulas. 

The M-dependence of integrated quantities like the energy for the PA and TIA, when 
M — > 00, is of the form 

1 \ 5 ( 1 x 5+2 



e = e ° + As {m) +Am {m) +•••• (20) 

with 5 = 2,4 for PA and TIA, respectively. The value 1/M in PIMC plays the same role 
than the step size h in a numerical integration and in both cases one is interested in the 
extrapolation to the ideal case h — (1/M — 0). Richardson extrapolation is a simple and 
clever way of performing that extrapolation improving the order of the approach. In the 
present case, the extrapolation to M = 00 is given by 

^-Mi^K?}**-*). (2l) 



Ei and E 2 being the energies estimated using Mi and M 2 beads, respectively (M 2 > Mi). 
When the calculation proceeds, and M is progressively increased it is also useful to know 
how far one is from the asymptotic law (20). In this respect, the Richardson extrapolation 
can also be used to predict the energy for a number of beads M > M 2 > Mi through 



III. RESULTS 

The usefulness of the Richardson extrapolation on top of the PIMC calculations and the 
accuracy of the TIA with respect to the PA have been studied in three different systems. The 
first one corresponds to the test problem of a particle in a one-dimensional harmonic potential 
in which the exact solution is known. In the second one, we study liquid Ne at 25.8 K, a 
real system at relatively high density and where the interatomic interaction presents a hard 
core at short distances. Finally, the most exigent test in the present analysis corresponds to 
the calculation of the energy of liquid 4 He at 5.1 K in which the low temperature, the small 
mass, and the hard core of the interactions make the quantum effects much larger. 

The sampling of the three systems studied has been carried out by combining collective 
movements of some beads of a given chain and movements of the center of mass of each 
one of the chains representing the atoms. In both movements the size of the proposed 
movements is fixed to keep a Metropolis acceptance of 30-50%. As it has been proved, the 
introduction of multibead movements is unavoidable when the number of beads increases 
since they eliminate the slowing down observed in a bead by bead sampling. To this end, 
we have used the bisection 2 ' 13 and staging 10-12 methods: the bisection up to level three for 
short chains and staging for longer ones. Both methods correspond to an exact sampling of 
the free (kinetic) part of the action and therefore, in the Metropolis step, only the potential 
part of the action is sampled. 

A. Harmonic oscillator 

We consider a particle in a one-dimensional harmonic well with Hamiltonian 




(22) 



H = 



2m dx 2 



H — moj x 
2 



(23) 
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In this problem, the partition function and the energy are exactly known, 6 

Z = [2sinh(/3/^/2)] _1 (24) 
E = -hLucoth{(3fiLj/2) . (25) 

In the PIMC calculation we have assumed u — h — m — 1, and T = 0.2. The exact energy 
from Eq. (25) at this temperature and using these reduced units is 0.50678. 

In Table I, PIMC results obtained using the PA and TIA approximations are reported 
as a function of the number of beads M. As already known for this model problem, 17 ' 18 
the accuracy of £tia is manifestly superior to the PA energies £pa- Within a statistical 
fluctuation of 10~ 5 , PA arrives to the exact energy when M = 512 whereas TIA requires 
only M = 32, a figure sixteen times smaller. In the same Table, we report Richardson 
extrapolations to a given M (E$ , Eq. 22) and to M — > oo (E^ , Eq. 21), derived from 
calculated results (-Epa, -Etia) using M/4 and M/2 beads. When 5 = 2, i.e., in the PA, the 
extrapolation to M — > oo achieves the asymptote with only M = 64, a factor of eight smaller 
than using only the direct output Epa- The extrapolation to a next M is also helpful to 



M 


£pa 


E m 

M 


E {2) 


Etia 


M 




2 


0.30755 






0.44702 






4 


0.43162 






0.50053 






8 


0.48424 


0.46264 


0.47298 


0.50630 


0.50387 


0.50410 


16 


0.50085 


0.49740 


0.50178 


0.50675 


0.50666 


0.50668 


32 


0.50528 


0.50500 


0.50639 


0.50678 


0.50678 


0.50678 


64 


0.50641 


0.50639 


0.50676 


0.50678 


0.50678 


0.50678 


128 


0.50669 


0.50669 


0.50679 




0.50678 


0.50678 


256 


0.50676 


0.50676 


0.50678 








512 


0.50678 


0.50678 


0.50678 








1024 




0.50678 


0.50679 









TABLE I: PA (-Epa) and TIA (.Etta) results for the one-dimensional harmonic oscillator at T = 0.2. 
Richardson extrapolations to M {E^) and to M — > oo (E^), using the M/4 and M/2 energies, 
are also reported (S = 2,4 stands for PA and TIA, respectively). 
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FIG. 1: PA (filled squares) and TIA (filled circles) results for the one-dimensional harmonic os- 
cillator as a function of 1/M. Solid and dashed lines are the fits (26) to the data. Open symbols 
stand for the Richardson extrapolations (squares, PA; circles, TIA). 

know when the expected quadratic law is attained since there that extrapolation coincides 
with the direct PA calculation. Finally, in Table I, the Richardson extrapolations from TIA 
results are shown. In this case, the gain achieved by the extrapolation is significantly smaller 
since only a factor of two seems reachable. 

The dependence of the PIMC results on the step 1/M is plotted in Fig. I. As it is 
expected from the accuracy of the PA and TIA factorizations, the behavior of the energy 
when 1/M — * follows the power law 



with E the exact energy and 5 = 2,4 for PA and TIA, respectively. The lines in the Figure 
correspond to fits to the PIMC results when the power law (26) starts to be valid. As 
commented before, Richardson extrapolation is a good tool to establish the smaller value of 
M from which up PIMC results follow this asymptotic behavior. On the other hand, one 
can see in the same Figure the behavior of the extrapolated energies. The combination of 
the PA and Richardson extrapolation reduces dramatically the difference between PA and 



E = E + As(l/M) 5 , 



(26) 
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M 


(V/N)p A 


(K/N) PA 


(E/N) PA 


(WAOtta 

V I /J- 1 r\ 


(K/N) T i A 


(£/iV) T TA 


2 


-210.25(12) 


44.42(2) 


-165.84(14) 


-207.18(13) 


48.04(4) 


-160.14(16) 


4 


-208.30(16) 


47.80(4) 


-160.50(19) 


-207.58(17) 


49.43(5) 


-158.15(22) 


8 


-207.80(17) 


49.05(5) 


-158.76(22) 


-207.48(11) 


49.61(3) 


-157.87(13) 


16 


-207.56(17) 


49.48(5) 


-158.08(20) 


-207.47(20) 


49.71(6) 


-157.76(22) 


32 


-207.52(22) 


49.60(6) 


-157.92(22) 








64 


-207.51(44) 


49.70(9) 


-157.81(36) 









TABLE II: Total (E), kinetic (K), and potential (V) energies for a different number of beads M 
in liquid Ne at 25.8 K. Figures within brackets are the statistical errors. 

TIA with zero computational cost. 
B. Neon at 25.8 K 

The first application to a real system corresponds to a PIMC simulation of liquid Ne at 
a experimental point of coordinates: temperature T = 25.8 K and density p = 0.0363 A~ 3 . 
We have used a simulation box containing 108 particles with periodic boundary conditions 
to represent the homogeneous liquid. Finite size effects have been studied and proved to 
be well corrected by summing up to the energy standard tail corrections to the potential 
energy. The interatomic potential between Ne atoms is the HFD-B model proposed by Aziz 
et al 30 

The results obtained using the PA and TIA for the total and partial energies are reported 
in Table II for an increasing number of beads. The value of the kinetic energy, which is 
compatible with previous estimations reported in Refs. 31 and 32, is a good measure of 
the relevance of quantum effects. Subtracting the classical kinetic energy (3T/2) from the 
PIMC value, this quantum correction amounts 11 K to be compared with the total value 
reported in Table II, 49.7 K . Taking into account the statistical error bars, the total energy 
reaches its asymptotic value in terms of 1/M for M = 32 and M = 8 using the PA and the 
TIA, respectively. Therefore, M is reduced by a factor of four, a figure significantly smaller 
than the gain obtained in the harmonic oscillator. In spite of this decrease, it is still more 
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M 


/ 7-1 / A7\ 


/ 7— I / i\t\ (2) 

( E / N ) M 


/ 7-1 / Ar^(2) 


/ 7-1 / A T \ 

(E/N) T i A 


/ 7-1 / 7V T \ (4) 

( e / n )m 


2 


-165.84(14) 






-160.14(16) 




4 


-160.50(19) 






-158.15(22) 




8 


-158.76(22) 


-159.16(24) 


-158.72(25) 


-157.87(13) 


-158.03(24) 


16 


-158.08(20) 


-158.32(28) 


-158.18(29) 


-157.76(22) 


-157.85(14) 


32 


-157.92(22) 


-157.91(25) 


-157.85(26) 




-157.75(24) 


64 


-157.81(36) 


-157.88(28) 


-157.87(29) 






128 




-157.78(45) 


-157.77(47) 







(E/N) 



(4) 



-158.02(24) 
-157.85(14) 
-157.75(24) 



TABLE III: Richardson extrapolations, to M ((E/N)$) and to oo ((E/N) ( £), of the total energy 
of liquid Ne with both, PA and TIA. 

efficient to use the TIA than the PA since the computation time for getting the same level 
of statistical error using TIA is less than the one required with PA. 

The usefulness of the Richardson extrapolation in this system is explored in Table III. 
Focusing the analysis on the PA action, in which the extrapolation is expected to be more 
useful, one can see how the extrapolated estimations approach the asymptote faster than the 
direct calculations: starting from M — 16 the energy remains constant inside the confidence 
interval determined by the statistical error. It is obvious that the unavoidable presence 
of statistical errors in a real calculation reduces the accuracy of the Richardson extrapola- 
tion. Nevertheless, it always improves the direct calculation and helps in determining the 
proximity to the expected power law (26). On the other hand, in the case of the TIA the 
extrapolation also improves the direct estimations but the gain is much smaller than in the 
PA. 

The dependence of the PA and TIA energies on the step 1/M is shown in Fig. 2. The lines 
correspond to fits (26) to the data when 1/M — > with 5 = 2,4 for PA and TIA, respectively. 
From the numerical fitting, the final energy per particle is -157.846(20) K for the PA and 
-157.824(42) K for the TIA. The Figure also contains the Richardson extrapolations to 
1/M — > from PA, which represent a significant improvement with respect the PA data and 
approach the fourth-order behavior of TIA. 
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-168 



-160 - 



-164 - 
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1/M 



0.4 
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FIG. 2: PA (filled squares) and TIA (filled circles) energies for Ne at 25.8 K as a function of 
1/M. Solid and dashed lines are the fits (26) to the data. Open squares stand for the Richardson 
extrapolations using the PA data. 

C. Liquid 4 He at 5.1 K 

The third system studied in the present work is the more demanding one: liquid 4 He 
at 5.1 K and density p = 0.02186 A~ 3 . As it is well known, liquid 4 He is the paradigm 
of a quantum Bose liquid and its properties have been well reproduced theoretically both 
at zero 33 and finite temperature. 2 At 5.1 K, 4 He is in its normal phase and the probability 
of accepting exchanges by introducing the correct symmetry in the action is very small, so 
its influence on the total energy becomes negligible. 2 Therefore, we consider 4 He atoms as 
Boltzmann-like particles. The interatomic potential is the HFD-B(HE) model from Aziz et 
a/., 34 used with high accuracy in zero-temperature calculations, and the homogeneous phase 
is simulated using a cubic box containing 64 atoms with periodic boundary conditions. 
Increasing the number of atoms, for a fixed number of beads, we have verified that the tail 
corrections added to the energy accounts satisfactorily for finite-size effects. 

Table IV contains results for the total and partial energies for an increasing number of 
beads M. The essentially quantum nature of 4 He is reflected in the high value of its kinetic 
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M 




(K/N)-px 


(E/N)pa 


tp A 


(WAOtta 


(K / AOtt a 


(E/N)tia 


6TT a 
v 1 1/A 


8 


-23.955(15) 


14.663(15) 


-9.292(20) 


0.8 


-22.803(10) 


16.503(25) 


-6.300(26) 


2.7 


16 


-22.737(10) 


16.885(20) 


-5.852(23) 


1.5 


-22.124(12) 


18.157(25) 


-3.967(26) 


4.5 


32 


-22.162(13) 


18.180(23) 


-3.982(27) 


3.3 


-21.895(12) 


18.852(26) 


-3.043(27) 


9.9 


64 


-21.940(16) 


18.754(23) 


-3.186(25) 


9.2 


-21.846(14) 


18.996(26) 


-2.850(27) 


26.1 


128 


-21.870(18) 


18.951(24) 


-2.919(26) 


17.1 


-21.850(14) 


19.027(27) 


-2.823(28) 


51.9 


256 


-21.848(17) 


19.006(23) 


-2.842(25) 


40.5 










512 


-21.843(20) 


19.036(30) 


-2.807(32) 


82.4 











TABLE IV: Total (E), kinetic (K), and potential (V) energies for a different number of beads 
M in liquid 4 He at 5.1 K. tpA (*tia) is the CPU time in thousands of seconds of a PIV-2.4GHz 
computer using the PA (TIA) approximation. 



M 


(E/N) PA 


{E/N)^ 


(E/N)W 


(E/N) T1A 




(E/N)^ 


8 


-9.292(20) 






-6.300(26) 






16 


-5.852(23) 






-3.967(26) 






32 


-3.982(27) 


-4.992(29) 


-4.705(31) 


-3.043(27) 


-3.821(28) 


-3.811(28) 


64 


-3.186(25) 


-3.514(34) 


-3.359(36) 


-2.850(27) 


-2.985(29) 


-2.981(29) 


128 


-2.919(26) 


-2.987(31) 


-2.921(33) 


-2.823(28) 


-2.838(29) 


-2.837(29) 


256 


-2.842(25) 


-2.852(33) 


-2.830(35) 




-2.821(30) 


-2.821(30) 


512 


-2.807(32) 


-2.823(31) 


-2.816(33) 








1024 




-2.798(40) 


-2.795(43) 









TABLE V: Richardson extrapolations, to M ((E/N) { ^) and to oo ((E/N) ( £>), of the total energy 
of liquid 4 He with both, PA and TIA. 

energy, comparable to the potential energy. The quantum correction on top of the classical 
value is ~ 11.5 K, to be compared with a total value of 19.0 K. Regarding the behavior of 
the PA and TIA energies with M, one can observe a similar behavior to the one obtained 
for Ne in the previous Subsection. The asymptotic value using PA is reached for M = 256 
whereas M = 64 for the TIA; the accuracy is then improved by a factor of four. Considering 
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-22.6 




17.0 1 1 1 1 1 

0.00 0.01 0.02 0.03 0.04 

1/M 



FIG. 3: Potential, total, and kinetic energies of liquid 4 He at 5.1 K. PA and TIA data are shown 
as filled squares and filled circles, respectively. Solid and dashed lines are the fits (26) to the data. 
Open squares stand for the Richardson extrapolations using the PA data. 

altogether, the reduction in the number of beads and the additional computation required 
for the TIA, the efficiency of the TIA is still greater than the one of the PA. 

Richardson extrapolations using the PIMC results for the total energies are reported in 
Table V. The main features are common to the ones observed in Ne. The use of Richardson 
extrapolation always improves the direct calculation with zero computational cost and it 
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arrives faster to the asymptote. From the results contained in the Table, one can see that 
the right energy using PA is obtained with data up to M = 128, a factor of two smaller than 
the direct PA calculation. In the case of the TIA, the extrapolation is also better than the 
direct output but it does not significantly improve the efficiency in achieving the plateau. 

Figure 3 is a plot of the total and partial energies as a function of 1/M near the region 
where the power laws (26) are satisfied. As it can be seen, the TIA data approaches the 
limit according to a fourth-order law, a point that sometimes has been questioned for hard- 
core-like interactions at very low temperatures. The Richardson extrapolations make a good 
job for both the total and partial energies by improving significantly the 1/M dependence of 
the PA data. An unbiased estimation of the energies is drawn from the numerical fits (26) 
with 5 = 2 for PA and 4 for TIA. Using PA, K/N = 19.026(7) K, V/N = -21.842(2) K, 
and E/N = -2.816(8) K. With the TIA, K/N = 19.019(6) K, V/N = -21.847(3) K, and 
E/N = —2.827(5) K. At the same thermodynamic point, Ceperley and Pollock 35 reported 
a total energy of -2.61 K and a kinetic energy of 18.09 K, but using a different interatomic 
potential. 

IV. CONCLUSIONS 

In the present work we have explored two possible strategies for improving the number-of- 
beads-dependence of the PIMC algorithm, specially for exigent problems in which the use of 
PA may require a very large number M. First, the usefulness of the Richardson extrapolation 
on top of the PA data has been analyzed for the first time. Richardson extrapolation is a 
clever way of practically improving the order of the errors introduced by discretization and 
its use in numerical algorithms like integration or solution of differential equations is widely 
known. Also in PIMC a discretization through an effective step 1/M is introduced and one 
is interested in going to the limit 1/M — > 0. The results obtained have proved that also in 
PIMC Richardson extrapolation can be useful since it approaches faster to the asymptote 
and helps to know when the discrete results behave according to the expected power law. It 
is worth mentioning that the presence of statistical errors mask somewhat the signal but it 
is also true that the computational cost of the extrapolation is certainly zero. 

The second main focus of the work has been to study the behavior of the TIA in systems 
with hard-core-like interactions and well into the quantum regime. Previous tests of the 
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TIA 24 ~ 27 were performed in less exigent conditions and some doubts concerning the fail of 
the predicted fourth-order accuracy were formulated. The present results, mainly the ones 
achieved in the calculation of 4 He, have shown that the fourth order law is maintained. 
On the other hand, the TIA introduces in the calculation additional computation due to 
the need of the second derivative of the potential (using the virial estimator for the kinetic 
energy, absolutely necessary when M increases) . This is not a serious problem, at least when 
empirical potentials are used, and the computer time to achieve the same statistical error 
in the asymptotic limit in M is ever reduced. Our present experience shows that the CPU 
time is at least a 30% smaller, and that this percentage increases with the number of beads 
required to reach the asymptote. 

As a final point, the fact that the double commutator [[V, K], V] in the TIA has allowed for 
a substantial decrease of the systematic errors for very different potentials is an encouraging 
result. At present, Chin 22 is developing various families of symplectic algorithms in which 
that double commutator originates a fourth-order correction which can be fine-tuned both 
in sign and magnitude. This bears the promise of achieving an effective cancellation of the 
leading fourth and sixth order corrections, feature which could result in a dramatic reduction 
of the computational cost of PIMC calculations. Further work in this direction would be 
most interesting. 
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